﻿///////////////////////////////////////////////////////////////////////////////
//
//  This file is part of MathLib.NET.
//
//  This library is free software; you can redistribute it and/or
//  modify it under the terms of the GNU Lesser General Public
//  License as published by the Free Software Foundation; either
//  version 2.1 of the License, or (at your option) any later version.
//  
//  This library is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//  Lesser General Public License for more details.
//  
//  You should have received a copy of the GNU Lesser General Public
//  License along with this library;  If not, see 
//  <http://www.gnu.org/licenses/>.
//
///////////////////////////////////////////////////////////////////////////////

using System;
using System.Collections.Generic;
using ILNumerics;

// Discrete Math functions
// =======================
//
// factor - Prime factors 
// factorial - Factorial function 
// gcd - Greatest common divisor 
// isprime - Array elements that are prime numbers 
// lcm - Least common multiple 
// nchoosek - Binomial coefficient or all combinations  
// perms - All possible permutations 
// primes - Generate list of prime numbers 
// rat, rats - Rational fraction approximation 


namespace MathLib
{
    public partial class MLMath
    {
        /// <summary>
        /// Object for performing integer factorizations.
        /// </summary>
        private static IntegerFactorize mIntFact = new IntegerFactorize();

        /// <summary>
        /// Factorial.
        /// </summary>
        /// <param name="n">Argument to return factorial of</param>
        /// <returns>Factorial of argument</returns>
        public static double factorial(int n)
        {
            double f;
            if (n >= 0)
            {
                double m = 1.0;
                f = 1.0;
                while (n != 0)
                {
                    f *= m;
                    m++;
                    n--;
                }
            }
            else
            {
                f = double.NaN;
            }

            return f;
        }

        /// <summary>
        /// Factorial.
        /// </summary>
        /// <param name="arg">Argument to return factorial of</param>
        /// <returns>Factorial of argument</returns>
        public static ILArray<double> factorial(ILArray<double> arg)
        {
            ILArray<double> facArray = new ILArray<double>(arg.Dimensions.ToIntArray());

            for (int i = 0; i < facArray.InternalArray4Experts.Length; i++)
            {
                facArray.InternalArray4Experts[i] = factorial((int)arg.InternalArray4Experts[i]);
            }

            return facArray;
        }

        /// <summary>
        /// Factorial.
        /// </summary>
        /// <param name="arg">Argument to return factorial of</param>
        /// <returns>Factorial of argument</returns>
        public static ILArray<double> factorial(ILArray<int> arg)
        {
            ILArray<double> facArray = new ILArray<double>(arg.Dimensions.ToIntArray());

            for (int i = 0; i < facArray.InternalArray4Experts.Length; i++)
            {
                facArray.InternalArray4Experts[i] = factorial(arg.InternalArray4Experts[i]);
            }

            return facArray;
        }

        /// <summary>
        /// Factorial.
        /// </summary>
        /// <param name="arg">Argument to return factorial of</param>
        /// <returns>Factorial of argument</returns>
        public static ILArray<double> factorial(ILArray<uint> arg)
        {
            ILArray<double> facArray = new ILArray<double>(arg.Dimensions.ToIntArray());

            for (int i = 0; i < facArray.InternalArray4Experts.Length; i++)
            {
                facArray.InternalArray4Experts[i] = factorial((int)arg.InternalArray4Experts[i]);
            }

            return facArray;
        }

        /// <summary>
        /// Find prime factors of integer.
        /// </summary>
        /// <param name="arg">Argument to factorize</param>
        /// <returns>Integer array containing result</returns>
        public static ILArray<int> factor(int arg)
        {
            List<int> pfactors = new List<int>();
            foreach (uint pfact in mIntFact.FindFactors((uint)arg))
            {
                pfactors.Add((int)pfact);
            }

            ILArray<int> factorArray = new ILArray<int>(new int[] { pfactors.Count });

            for (int i = 0; i < pfactors.Count; i++)
            {
                factorArray.InternalArray4Experts[i] = pfactors[i];
            }

            mIntFact.ClearPrimes();

            return factorArray;
        }

        /// <summary>
        /// Find prime factors of integer.
        /// </summary>
        /// <param name="arg">Argument to factorize</param>
        /// <returns>Integer array containing result</returns>
        public static ILArray<uint> factor(uint arg)
        {
            List<uint> pfactors = new List<uint>();
            foreach (uint pfact in mIntFact.FindFactors(arg))
            {
                pfactors.Add(pfact);
            }

            ILArray<uint> factorArray = new ILArray<uint>(new int[] { pfactors.Count });

            for (int i = 0; i < pfactors.Count; i++)
            {
                factorArray.InternalArray4Experts[i] = pfactors[i];
            }

            mIntFact.ClearPrimes();

            return factorArray;
        }

        /// <summary>
        /// Greatest common divisor of each element of matrix a and b.
        /// </summary>
        /// <param name="mata">First matrix</param>
        /// <param name="matb">Second matrix</param>
        /// <returns>Greatest common divisor</returns>
        public static ILArray<int> gcd(ILArray<int> mata, ILArray<int> matb)
        {
            ILArray<int> gcdArray = new ILArray<int>(mata.Dimensions.ToIntArray());

            for (int i = 0; i < gcdArray.InternalArray4Experts.Length; i++)
            {
                gcdArray.InternalArray4Experts[i] = (int)MathNet.Numerics.NumberTheory.IntegerTheory.GreatestCommonDivisor((long)mata.InternalArray4Experts[i], (long)matb.InternalArray4Experts[i]);
            }

            return gcdArray;
        }

        /// <summary>
        /// Greatest common divisor of each element of matrix a and b.
        /// </summary>
        /// <param name="mata">First matrix</param>
        /// <param name="matb">Second matrix</param>
        /// <returns>Greatest common divisor</returns>
        public static ILArray<uint> gcd(ILArray<uint> mata, ILArray<uint> matb)
        {
            ILArray<uint> gcdArray = new ILArray<uint>(mata.Dimensions.ToIntArray());

            for (int i = 0; i < gcdArray.InternalArray4Experts.Length; i++)
            {
                gcdArray.InternalArray4Experts[i] = (uint)MathNet.Numerics.NumberTheory.IntegerTheory.GreatestCommonDivisor((long)mata.InternalArray4Experts[i], (long)matb.InternalArray4Experts[i]);
            }

            return gcdArray;
        }

        /// <summary>
        /// Lowest common multiple of each element of matrix a and b.
        /// </summary>
        /// <param name="mata">First matrix</param>
        /// <param name="matb">Second matrix</param>
        /// <returns>Lowest common multiple</returns>
        public static ILArray<int> lcm(ILArray<int> mata, ILArray<int> matb)
        {
            ILArray<int> lcmArray = new ILArray<int>(mata.Dimensions.ToIntArray());

            for (int i = 0; i < lcmArray.InternalArray4Experts.Length; i++)
            {
                lcmArray.InternalArray4Experts[i] = (int)MathNet.Numerics.NumberTheory.IntegerTheory.LeastCommonMultiple((long)mata.InternalArray4Experts[i], (long)matb.InternalArray4Experts[i]);
            }

            return lcmArray;
        }

        /// <summary>
        /// Lowest common multiple of each element of matrix a and b.
        /// </summary>
        /// <param name="mata">First matrix</param>
        /// <param name="matb">Second matrix</param>
        /// <returns>Lowest common multiple</returns>
        public static ILArray<uint> lcm(ILArray<uint> mata, ILArray<uint> matb)
        {
            ILArray<uint> lcmArray = new ILArray<uint>(mata.Dimensions.ToIntArray());

            for (int i = 0; i < lcmArray.InternalArray4Experts.Length; i++)
            {
                lcmArray.InternalArray4Experts[i] = (uint)MathNet.Numerics.NumberTheory.IntegerTheory.LeastCommonMultiple((long)mata.InternalArray4Experts[i], (long)matb.InternalArray4Experts[i]);
            }

            return lcmArray;
        }

        /// <summary>
        /// Determines if each matrix element is prime.
        /// </summary>
        /// <param name="mat">Matrix</param>
        /// <returns>1 for prime elements, 0 for non-prime elements</returns>
        public static ILArray<int> isprime(ILArray<int> mat)
        {
            ILArray<int> ipArray = new ILArray<int>(mat.Dimensions.ToIntArray());

            for (int i = 0; i < ipArray.InternalArray4Experts.Length; i++)
            {
                ipArray.InternalArray4Experts[i] = mIntFact.IsPrime((uint)mat.InternalArray4Experts[i]) ? 1 : 0;
            }

            mIntFact.ClearPrimes();

            return ipArray;
        }

        /// <summary>
        /// Determines if each matrix element is prime.
        /// </summary>
        /// <param name="mat">Matrix</param>
        /// <returns>1 for prime elements, 0 for non-prime elements</returns>
        public static ILArray<uint> isprime(ILArray<uint> mat)
        {
            ILArray<uint> ipArray = new ILArray<uint>(mat.Dimensions.ToIntArray());

            for (int i = 0; i < ipArray.InternalArray4Experts.Length; i++)
            {
                ipArray.InternalArray4Experts[i] = (uint)(mIntFact.IsPrime(mat.InternalArray4Experts[i]) ? 1 : 0);
            }

            mIntFact.ClearPrimes();

            return ipArray;
        }

        /// <summary>
        /// Returns the number of rows produced for given values of n and k.
        /// </summary>
        /// <param name="n">Number of elements in row vector</param>
        /// <param name="k">Number of items taken at a time</param>
        /// <returns>Number of combinations</returns>
        public static double nchoosek(int n, int k)
        {
            return factorial(n) / (factorial(n - k) * factorial(k));
        }

        /// <summary>
        /// Returns the number of rows produced for given values of n and k.
        /// </summary>
        /// <param name="n">Number of elements in row vector</param>
        /// <param name="k">Number of items taken at a time</param>
        /// <returns>Number of combinations</returns>
        public static ILArray<int> nchoosek(ILArray<int> n, int k)
        {
            ILArray<int> nckArray = new ILArray<int>(new int[] { k, (int)nchoosek(n.InternalArray4Experts.Length, k) });

            int[] idx = new int[k];
            for (int i = 0; i < k; i++) idx[i] = i;

            int row = 0;
            while (true)
            {
                // select items
                for (int i = 0; i < k; i++) nckArray.InternalArray4Experts[row * k + i] = n.InternalArray4Experts[idx[i]];

                // try to increment index
                bool incremented = false;
                int m = n.InternalArray4Experts.Length - 1;
                for (int i = 0; i < k; i++)
                {
                    // if there is space to move right
                    if (idx[k - i - 1] < m)
                    {
                        // increment index and stop searching
                        int j = k - i - 1;
                        idx[j++]++;
                        for (; j < k; j++) idx[j] = idx[j - 1] + 1;
                        incremented = true;
                        row++;
                        break;
                    }
                    m--;
                }

                // if no increment was possible then all combinations have been found
                if (!incremented)
                    break;
            }

            return nckArray;
        }

        /// <summary>
        /// Generate all permutations from a row vector.
        /// </summary>
        /// <param name="array">Array to generate permutations for</param>
        /// <returns>Array containing permutations of original vector</returns>
        public static ILArray<int> perms(ILArray<int> array)
        {
            ILArray<int> permArray = new ILArray<int>(new int[] { array.InternalArray4Experts.Length, (int)factorial(array.InternalArray4Experts.Length) });

            perms(array, permArray, array.InternalArray4Experts.Length);

            return permArray;
        }

        public static void perms(ILArray<int> array, ILArray<int> permArray, int n)
        {
            if (n > 1)
            {
                perms(array, permArray, n - 1);

                // add new element
                int cols = array.InternalArray4Experts.Length;
                for (int i = 0; i < n - 1; i++)
                {
                    // p[n - 1][i] = a[n - 1]
                    permArray.InternalArray4Experts[i * cols + n - 1] = array.InternalArray4Experts[n - 1];
                }

                // generate permutations
                for (int i = 1; i < n; i++)
                {
                    for (int j = 0; j < n; j++)
                    {
                        for (int k = 0; k < n - 1; k++)
                        {
                            // p[(j + 1) % n][i * (n - 1) + k] = p[j][(i - 1) * (n - 1) + k]
                            permArray.InternalArray4Experts[(i * (n - 1) + k) * cols + (j + 1) % n] = permArray.InternalArray4Experts[((i - 1) * (n - 1) + k) * cols + j];
                        }
                    }
                }
            }
            else
            {
                permArray.InternalArray4Experts[0] = array.InternalArray4Experts[0];
            }
        }

        /// <summary>
        /// Returns a row-vector of prime numbers less than or equal to n
        /// </summary>
        /// <param name="n">Maximum prime number</param>
        /// <returns>Vector of primes up to maximum inclusive</returns>
        public static ILArray<uint> primes(int n)
        {
            foreach (uint prime in mIntFact)
                if (prime > n) break;

            IList<uint> primeList = mIntFact.GetPrimeList();
            ILArray<uint> primArray = new ILArray<uint>(new int[] { 1, primeList.Count });

            for (int i = 0; i < primeList.Count; i++)
                primArray.InternalArray4Experts[i] = primeList[i];

            mIntFact.ClearPrimes();

            return primArray;
        }

        /// <summary>
        /// Returns if the rational number represented by the specified d-coefficients is
        /// accurate to the specified tolerance.
        /// </summary>
        /// <param name="d">D-coefficients</param>
        /// <param name="target">Target number</param>
        /// <param name="tolerance">Tolerance</param>
        /// <returns>true if accurate within tolerance, otherwise false</returns>
        private static bool ratisaccurate(IList<int> d, double target, double tolerance)
        {
            if (d.Count < 1) return target == 0.0;

            double val = d[d.Count - 1];
            for (int i = 1; i < d.Count; i++)
            {
                val = 1 / val;
                val += d[d.Count - i - 1];
            }

            if (Math.Abs(target - val) < tolerance)
                return true;

            return false;
        }

        /// <summary>
        /// Gets the string that represents the given continued fraction coefficients.
        /// </summary>
        /// <param name="d">Continued fraction coefficients</param>
        /// <returns>String representation of fraction</returns>
        private static string ratfracstr(IList<int> d)
        {
            if (d.Count < 1) return "0";

            string frac = d[d.Count - 1].ToString();
            for (int i = 1; i < d.Count; i++)
            {
                string subfrac = i == 1 ? frac : "(" + frac + ")";
                frac = d[d.Count - i - 1].ToString() + " + 1 / " + subfrac;
            }

            return frac;
        }

        /// <summary>
        /// Calculates D-coefficients for the fractional approximation to the
        /// given argument at the specified accuracy.
        /// </summary>
        /// <param name="arg">Floating point value to approximate</param>
        /// <param name="tolerance">Tolerance of approximation accuracy</param>
        /// <returns>D coefficients representing fraction</returns>
        private static IList<int> CalcDCoeffs(double arg, double tolerance)
        {
            List<int> d = new List<int>();

            double val = arg;
            while (!ratisaccurate(d, arg, 1E-6))
            {
                d.Add((int)Math.Floor(val));

                val = 1.0 / (val - Math.Floor(val));
            }

            if (d[d.Count - 1] == 1 && d.Count > 1)
            {
                d.RemoveAt(d.Count - 1);
                d[d.Count - 1]++;
            }

            return d;
        }

        /// <summary>
        /// Gets the numerator and denominator of a fraction that
        /// approximates the value of arg to the specified tolerance.
        /// </summary>
        /// <param name="n">Numerator</param>
        /// <param name="k">Denominator</param>
        /// <param name="d">D-coefficients</param>
        /// <param name="tolerance">Minimum tolerance</param>
        private static void simplifycfrac(out int n, out int k, IList<int> d, double tolerance = 1E6)
        {
            if (d.Count == 0)
            {
                n = 0;
                k = 1;
                return;
            }

            n = 1;
            k = d[d.Count - 1];

            for (int i = 1; i < d.Count; i++)
            {
                // Add fractions
                n += d[d.Count - i - 1] * k;

                // Swap n and k if not last
                if (i < d.Count - 1)
                {
                    int oldk = k;
                    k = n;
                    n = oldk;
                }
            }
        }

        /// <summary>
        /// Output the continued fraction for the specified argument.
        /// </summary>
        /// <param name="arg">Floating point number to approximate</param>
        public static void rat(double arg)
        {
            IList<int> d = CalcDCoeffs(arg, 1E6);

            System.Diagnostics.Debug.WriteLine(ratfracstr(d));
        }

        /// <summary>
        /// Generates numerator and denominator of fractional approximation to arg.
        /// </summary>
        /// <param name="n">Numerator</param>
        /// <param name="k">Denominator</param>
        /// <param name="arg">Floating point number to approximate</param>
        /// <param name="tolerance">Minimum tolerance of approximation</param>
        public static void rat(out ILArray<int> n, out ILArray<int> k, ILArray<double> arg, double tolerance)
        {
            n = new ILArray<int>(arg.Dimensions.ToIntArray());
            k = new ILArray<int>(arg.Dimensions.ToIntArray());

            for (int i = 0; i < arg.InternalArray4Experts.Length; i++)
            {
                IList<int> d = CalcDCoeffs(arg.InternalArray4Experts[i], tolerance);

                simplifycfrac(out n.InternalArray4Experts[i], out k.InternalArray4Experts[i], d, tolerance);
            }
        }

        /// <summary>
        /// Generates numerator and denominator of fractional approximation to arg.
        /// </summary>
        /// <param name="n">Numerator</param>
        /// <param name="k">Denominator</param>
        /// <param name="arg">Floating point number to approximate</param>
        /// <param name="tolerance">Minimum tolerance of approximation</param>
        public static void rat(out ILArray<int> n, out ILArray<int> k, ILArray<double> arg)
        {
            rat(out n, out k, arg, 1E6);
        }

        /// <summary>
        /// Generates a string containing the fractional approximation to the supplied floating point value.
        /// </summary>
        /// <param name="arg">Value to approximate</param>
        /// <returns>Fractional approximation to floating point value</returns>
        public static string rats(double arg)
        {
            IList<int> d = CalcDCoeffs(arg, 1E6);

            int n, k;
            simplifycfrac(out n, out k, d);

            return string.Format("{0} / {1}", n, k);
        }
    }
}
